The Flint Hills of Kansas and Oklahoma are the largest remaining tracts of tallgrass prairie in North America. This area has undergone major changes in land management practices in the past 30 years. Traditional season-long cattle stocking with variable burn schedules has diversified to include intensive-early cattle stocking accompanied by annual burning. To understand how different land management practices affect the herpetofauna of a tallgrass prairie, we used mark-recapture statistics to analyze herpetofaunal community dynamics. We analyzed survey data collected over a 15-year time span (1989–2003) from a rangeland site in Cowley County, KS, USA. A modified Jolly-Seber open population model, POPAN-5, was used to estimate four community parameters: probability of species loss (φ;′), probability of detection (P), probability of entry (Pent), and species richness (N). The top models included burn status as a covariate for species loss rate, while cattle stocking received moderate support as a covariate. Rates of species loss were higher during burn years (φ;′ = 0.04, 95% CI: 0.02 to 0.08) than nonburn years (φ;′ = 0.00, 95% CI: 0.00 to 0.01). Analysis of the impacts of different management practices was difficult due to confounding effects of changes in both burning and grazing. Declines in species richness tended to be steepest during a period of season-long stocking, but results were not statistically significant. Though our limited data set does not allow us to draw strong conclusions on the effects of land management on herpetofaunal populations, the mark-recapture models illustrated in our study should prove to be a valuable tool in future analyses of similar data.
The tallgrass prairie has been reduced to less than five percent of its pre-European settlement extent, making it the most heavily impacted ecosystem in North America (Samson and Knopf, 1994). The largest contiguous section of remaining tallgrass prairie occurs in the Flint Hills of northern Oklahoma and Kansas (Knapp and Seastedt, 1998). Much of this region is unsuitable for row-crop agriculture because shallow soils are associated with limestone bedrock. However, the Flint Hills provide extensive rangeland for cattle. Much of the rangeland experienced a marked change in land management in the early 1980s. Traditionally, ranchers practiced season-long cattle stocking, which allowed a limited number of cattle to graze year-round, accompanied by rotation of burned and unburned areas every 2–3 years. Under intensive-early stocking, cattle stocking rates were roughly double that of the previous practice for only 90–120 days in the growing season with annual spring burns, which increased daily steer gains and gains per acre (Robbins et al., 2002; Smith and Owensby, 1978). Intense cattle grazing has been found to have a significant effect on the diversity of native plant species, reducing overall cover of indiangrass (Sorghastrum nutans), and increasing Kentucky bluegrass (Poa pratensis) and big bluestem (Andropogon gerardii) basal area (Owensby et al., 1988).
Fire is an essential component in the development and persistence of the tallgrass prairie ecosystem (Axelrod, 1985). Regular burning of the tallgrass prairie increases plant productivity, decreases aboveground biomass by removing dead litter, and limits woody vegetation encroachment (Briggs and Knapp, 1995; Heisler et al., 2003; Hulbert, 1988). However, annual burning eliminates aboveground cover, which may increase short-term predation risk for larger-bodied species of snakes (Wilgers, 2005).
Species richness estimates can be used as an indicator of disturbance effects on communities by comparing the numbers of species pre- and post-disturbance within in a single area. Past studies have determined that reptile and amphibian abundance and community composition are influenced by fire in grassland and shrub ecosystems, including desert (Fyfe, 1980), sandhill (Mushinsky, 1985), sand-pine scrub (Greenberg et al., 1994), spinifex grassland (Masters, 1996; Pianka, 1996), and tropical savanna (Braithwaite, 1987; Griffiths and Christian, 1996; Hannah and Smith, 1995; Trainor and Woinarski, 1994), as well as other ecosystems in North America (Pilliod et al., 2003; Russell et al., 1999). In the tallgrass prairie, past projects have studied the effects of fire on arthropods (Evans, 1984, 1988; Joern, 2004, 2005; Seastedt, 1984a,b), birds (Beutler, 1993; Powell, 2004; Zimmerman, 1992, 1997), and mammals (Clark and Kaufman, 1990; Clark et al., 1989; Kaufman et al., 1990; McMillan et al., 1995). A limited number of studies have examined the distribution of herpetofauna in tallgrass prairie (Busby and Parmelee, 1996; Busby et al., 1994; Heinrich and Kaufman, 1985), but the effects of fire and grazing on herpetofauna inhabiting the tallgrass prairie ecosystem are poorly understood (Cavitt, 2000; Setser and Cavitt, 2003).
To estimate species richness of an area (N ˆ), most studies have used unadjusted count data (Braithwaite, 1987; Fredericksen and Fredericksen, 2002; Hannah and Smith, 1995), Shannon-Wiener indices (Chao and Shen, 2003; Greenberg et al., 1994), or species accumulation curves (Thompson et al., 2003). All of these techniques assume that capture probabilities (P) equal one. Previous mark-recapture studies of amphibians and reptiles show that P < 1 for both individuals and cryptic species (Bailey et al., 2004a,b,c; Schmidt, 2003; Schmidt et al., 2002), which can result in a biased estimator of N ˆ (Pollock et al., 1990). Mark-recapture statistics can be used to analyze a variety of data types, which include but are not restricted to the typical mark-recapture methods, which utilize the marking of individuals. Mark-recapture models have been used to estimate species richness in spiders (Follner and Henle, 2001) and birds (Boulinier et al., 1998; Nichols et al., 1998), but not reptile or amphibian species. This is the first study to examine community dynamics of reptiles and amphibians using mark-recapture models. Our goals in this study were two-fold: (1) to estimate species richness and probabilities of detection for reptiles and amphibians, and (2) to determine the effects of land management techniques on community dynamics.
Materials and Methods
Study Site
Fourteen herpetofaunal surveys were conducted during a 15-year period (1989–1990 and 1992–2003) at Bud Jan Nitschke (BJN) cattle ranch in Cowley County, Kansas (37°15′ N, 96°43′ W). The 650 ha ranch is located in the southern portion of the Flint Hills, and the vegetation is dominated by a matrix of perennial warm-season C4 grasses, such as big bluestem (A. gerardii), little bluestem (A. scoparius), indiangrass (S. nutans), and switchgrass (Panicum virgatum). A diverse mixture of other less abundant species includes warm-season and cool-season grasses, composites, legumes, and other forbs. Woody species such as leadplant (Amorpha canescens), buckbrush (Symphoricarpos orbiculatus), New Jersey tea (Caenothus herbaceous), and smooth sumac (Rhus glabra) are locally common (Freeman and Hulbert, 1985).
Average monthly temperature ranges from a low in January (−2.7 C) to a high in July (26.6 C). Average annual total precipitation ranges from ~750 to 850 mm with 75% falling during a summer growing season (March–July; Bark, 1987). One common characteristic of the region is protruding limestone outcrops on hilltops, which provide suitable habitat for several species of amphibians and reptiles.
During the 15-year study, two management practices were implemented at the BJN ranch. Traditional season-long stocking (200 hundred cows with their calves, ~0.6 animals/ha) was implemented the first 10 years (1989–1998) and pastures were burned in alternate years. Starting in 1999, the land management practice switched to intensive-early cattle stocking (approximately 650 yearlings, ~1.0 animals/ha, maintained for 3 months starting in late spring) combined with annual burning.
Field Protocol
Herpetofaunal surveys were conducted on one day between 19–30 April each year. The survey route consisted of the same 4-km walk which took approximately four hours (1000–1400 h). Walking surveys is an effective and widespread practice for monitoring amphibian and reptile communities (Doan, 2003; Woodford and Meyer, 2003). The route included normal tallgrass prairie habitats, the shorelines of a pond, and intermittent streams. Amphibians and reptiles were detected by turning over rocks and other debris, and by sighting animals in the open. Each year, we recorded the number of individuals of each species encountered, and environmental conditions including: midday air temperature (C), water temperature (C), cloud cover, burn conditions (burned versus unburned), and the presence of flowing water in the intermittent streams (Table 1).
Table 1
Burn status and climatic conditions of surveys at Cowley County site.
Statistical Analysis
Mark-recapture models can be used to analyze a wide variety of data, which include but are not restricted to the more traditional mark-recapture methods requiring the marking of single individuals. More recent mark-recapture statistics enable us to estimate community parameters such as species richness, emigration, and immigration via species presence/absence data across repeated sampling efforts with no individual marking. The reference of mark-recapture statistics and models throughout our paper does not imply the use of traditional mark-recapture methods.
Species richness is used to compare the complexity of communities between two or more habitats or to compare before and after a perturbation. When measuring community dynamics, two additional parameters are of interest: (1) the probability of species loss (local extinction or permanent emigration), and (2) the probability of a new species entering the community (immigration). Mark-recapture statistics were used to estimate these initial parameters along with the species richness of the herpetofaunal community of the study site. Initial detection data were transformed into presence/absence data per species by year, creating a distinct encounter history for each of the species recorded across all 14 years of data (Table 2).
Table 2
Species detected during herpetofaunal surveys at Cowley County site (1989–2003).
Program POPAN-5, a Cormack-Jolly Seber (CJS) model design was fit to the data structure (Schwarz and Arnason, 1996). This model assumes an open population (allows immigration, emigration, births, and deaths in the community) throughout the entire study period with near instantaneous sampling periods. Four parameters are estimated by POPAN-5: φ;, species retention; P, probability of any species detection given the species is in the area; Pent, probability of a new species entering the community; and N, species richness. Rates of species loss (φ;′) were calculated as the complement of retention (1 − φ;). Two other parameters can be derived from the parameter estimates above: b, number of species entries into the community (Equation 1); and N ˆ, species richness estimates for each individual time step (i) within the study starting at time zero with N = 0 (Equation 2).
Mark-recapture analyses were conducted using Program MARK (Version 4.1; White and Burnham, 1999) in three-steps. First, we selected a global model. Species retention (φ;), probability of detection (P), and probability of entry (Pent) were all modeled with time dependency (φ;t, Pt, Pent t, N) because we expected variation in annual conditions might impact these parameters. This was the most parameterized model (no. of parameters = 43).
Second, we tested the global model for possible overdispersion. Program RELEASE was used to estimate a variance inflation factor (ĉ), which was calculated by dividing the sum of all the chi-square tests by the sum of the degrees of freedom over all the component tests of Test 2 and 3 (Burnham et al., 1987).
Last, we fit reduced models where annual variation in the parameter estimates was modeled as a function of environmental covariates or as a constant. We modeled φ; as a function of burn status because spring fires are known to cause direct mortality of reptiles (Erwin and Stasiak, 1979; McLean, 1990). Probability of detection (P) was modeled as a function of temperature and cloud cover because amphibians and reptiles are ectothermic and may have reduced activity or seek refugia under cool conditions. Burn status and cloud cover were treated as Bernoulli conditions (1 = yes, 0 = no). In POPAN-5, species richness (N) is estimated as a single parameter, and the derived N per occasion cannot be modeled as a function of a time dependent environmental covariate. In surveys of vertebrate communities, most species are detected in the first few years, and new detections become progressively more difficult as time goes by. Thus, the probability of entry of new species into the community is likely to decline, and we modeled Pent as a linear function.
Models were constructed using the following link functions: logit (φ; and P), Mlogit1 (Pent), and log (N). The logit link constrains species retention and probability of detection from zero to one, which is necessary because both parameters are probabilities. The Mlogit1 link constrains the sum of all time steps (i.e., 1–14) from zero to one, which is necessary for a probability of entry (Pent) estimate due to its similarities with detection functions. The log link function does not constrain the estimate within a bounded, which is desired because estimate of species richness (N) are >1.
All possible combinations between our global model and the environmental covariates yielded a candidate set of 25 models. Model fit was assessed using quasi-Akaike's Information Criteria (QAICc), which corrects for overdispersion and small sample size (Burnham and Anderson, 1998). Model selection was based on the difference between QAICc values (ΔQAICc). The best-fit model had a ΔQAICc of zero, but models with a ΔQAICc ≤ 2 were considered equally parsimonious.
Parameter estimates, their corresponding standard errors, and 95% confidence intervals for the top models were calculated in Program MARK using all models. Derived estimates of N for each year were averaged across all models using the model averaging procedure of Program MARK. In a post-hoc analysis, we looked at the effect of land management practices on species richness by pooling years to give species richness estimates for three blocks of time: (I) initial years (1989–1992), (II) season-long stocking (1993–1998), and (III) intensive-early stocking (1999–2003). We did not include the initial years (1989–1992) in the land management estimates because high probability of species entry during the first three surveys indicated that initial surveys missed >5% of the species in the area. This model was not used for model selection or parameter estimates. Significant differences in species loss rates were tested using program CONTRAST (Hines and Sauer, 1989). Significance testing was performed using t-tests assuming unequal variances, and α was set at 0.05.
Results
Over the 14 years of this study, 1563 individuals representing 35 species of amphibians and reptiles were detected at the study site. The sample included species representing 4 orders of heterothermic tetrapods (Caudata, Anura, Testudinata, and Squamata) with two suborders of squamates (Sauria and Serpentes; Table 2). Components of Program RELEASE (Test 2 + Test 3) yielded a moderate variance inflation factor (χ2 = 64.7, df = 23, ĉ = 2.81) for the global model (φ;t, Pt, Pent t, N), indicating minor overdispersion of the count data.
Model Selection
The best fit models included constant parameters and the environmental covariate, burn status. The two top models had good support (wi > 0.20) and were considered to be equally parsimonious (ΔQAICc ≤ 2), while eight other models containing the additional environmental covariates (cattle stocking, air temperature, and cloud cover; Table 1) showed a moderate fit to the data (ΔQAICc ≤ 7, wi = 0.01 to 0.10; Table 3). The global model was a poor fit to the data (φ;t, Pt, Pent t, N; ΔQAICc = 79.9). The 3 time block model (φ;3tb, Pc, Pent, N), model 3, was only moderately supported by the data (ΔQAICc = 4.16, wi = 0.033).
Table 3
Summary of model selection for species richness estimates of reptiles and amphibians in tallgrass prairie, Cowley County, Kansas.
Of the candidate models, φ; was best modeled as a constant, or with the burn status covariate. Land management showed moderate support with ΔQAICc = 2.1 and garnished 10% of the data's support (Table 3). Probability of detection (P) was best modeled as a constant or with the two environmental covariates: mid-day air temperature and cloud cover (Table 3). Modeling Pent as a linear function was a better fit to the data than time dependence (ΔQAICc ≥ 20; Table 3).
Community Dynamics
Community demographics φ; and P were taken from the top two equally parsimonious models. The φ; parameter for model 1 yielded 1 rate across all 14 years (φ;′ = 0.025, 95% CI = 0.011 to 0.048), while model 2 yielded 2 rates according to burn status. Species loss rates were significantly higher during burn years (φ;′ = 0.04, 95% CI: 0.02 to 0.08) when compared to years without burning (φ;′ = 0.0, 95% CI = 0.00 to 0.01; Χ;2 = 8.16, df = 1, P = 0.004). Probability of detection was moderate in the top two models: model 1, P = 0.581 (95% CI = 0.524 to 0.636), and model 2, P = 0.586 (95% CI = 0.529 to 0.640). Probability of species entry into the community was best modeled as a linear function over time and declined from Pent = 0.19 in 1989 to Pent ≤ 0.03 by 1993 (Fig. 1a).
Figure 1
(A) Probability of new herpetofaunal species entering Cowley County tallgrass prairie site from 1989–2002.
Parameter estimates (x¯ ± 1 SE) were calculated by program MARK through model averaging of all candidate models (Table 1). (B) Estimates of all herpetofaunal species richness from 1989–2003 at Cowley County tallgrass prairie site. For post-hoc analysis, the surveys were divided into three time periods: (I) Initial years (1989–1992) where Pent ≥ 0.05, (II) Traditional season-long stocking (1993–1998), (III) Intensive-early stocking (1999–2003). Parameter estimates (x¯ ± 1 SE) were calculated by three time blocks model. Overall estimate (N ˆ) was calculated across all 14 years of study.
The overall estimate for species richness (N ˆ), across the entire time series did not significantly differ between the two top models, model 1 (N ˆ = 35.1) and model 2 (N ˆ = 35.0). These estimates matched the total number of species observed throughout the study. Annual estimates of species richness are a derived parameter and were averaged across all candidate models and years to estimate the average number of species in the community each year. Species richness estimates for each individual year did not reach the overall estimate across all years. Yearly species richness estimates were at a maximum in 1995 (N ˆ = 32.0, 95% CI = 29.4 to 34.5), with subsequent declines by the end of the study in 2003 (N ˆ = 26.6, 95% CI = 21.0 to 32.1; Fig 1b).
In the three block model, estimates of species loss rates were not different between the three time periods (φ;′I = 0.02, 95% CI = 0.0 to 0.08; φ;′II = 0.03, 95% CI = 0.01 to 0.08; φ;′III = 0.01, 95% CI = 0.01 to 0.39). Species richness estimates declined faster during the season-long stocking compared to the intensive-early stocking period (Fig. 1b), however, the variability in species richness estimates are high and greatly overlap.
Discussion
This is the first study to use mark-recapture statistics to analyze effects of land management practices on the community dynamics of herpetofauna. Overall, detection rates of herpetofaunal species were low throughout the study (P < 0.59). Mark-recapture statistics were a useful approach here because not all herpetofaunal species were detected every year and the probability of detection was affected by environmental covariates (climatic conditions), (Schmidt et al., 2002). Analyses assuming probability of detection close to unity (P = 1) would have resulted in biased estimates of N ˆ (Schmidt, 2003). Some of the annual variation in detection rates was due to environmental variables such as air temperature and cloud cover.
In this study, total species richness of the herpetofaunal community was estimated at 35 species across a 14-year time series in southeast Kansas. By comparison, 57 herpetofaunal species have been reported from Cowley County, and 97 species are on the current state list for Kansas (Collins, 1993; Potts and Collins, 2005). According to the probability of species entry into the system (Fig. 1a), our sampling detected a majority of species occurring in the area within the first three years of the study.
Effects of Management Practices
This retrospective analysis limits our ability to draw strong conclusions about the effects of land management on prairie herpetofauna. One sampling effort per year may not sufficiently explain species presence or absence in the area due to variations in activity based on climate or other variables. The trends of herpetofaunal community dynamics found in this study may be widespread throughout the area, but the lack of a reference sample does not allow us to determine this.
However, based on our limited data, the mark-recapture models we used suggest that different land management practices may have an effect on herpetofaunal species richness in the Flint Hills ecosystem. Intensive-early stocking, which is a combination of repeated annual burning and intense grazing, has been shown to have a negative effect on other tallgrass prairie fauna, such as birds (Robbins et al., 2002; Rohrbaugh et al., 1999; Zimmerman, 1997) and mammals (Kaufman and Kaufman, 1997).
Species loss rates were significantly higher in years after a spring burn. During the first 10 years of season-long stocking, species richness peaked around 31. Species loss rates declined during both management practices. Declines were steepest during season-long stocking, however variability in species richness estimates were greater during intensive-early stocking. The management practice and three time block models were only given moderate support from the data, and were not among the top models. Failure of the three time block model to improve the model fit indicates minimal influence of the initial detection time period on the overall community dynamics estimates. Possible confounding effects between the observed negative impact of burning and the unknown impact of grazing may have decreased the overall impact of the two stocking practices, thus weakening the model fit to the data. Thus, management practice may influence herpetofaunal community dynamics, but, further research is necessary to determine any significant differences.
Herpetofaunal species loss rates were significantly higher following burn years across the entire 14-year data set. Fire and grazing disturb the vegetation structure similarly, leaving much of the disturbed area devoid of grass more than a few centimeters high for much of the spring (Robbins et al., 2002). Intensive-early stocking provides grasses a chance to recover during later months of the year. This results in a continuous fuel source and more intense fires (Smith and Owensby, 1978). Direct mortality of herpetofauna due to fire can have large effects on certain species within the herpetofaunal community (Erwin and Stasiak, 1979; McLean, 1990), and more intense fires can lead to an increase in direct mortality of individuals in the area. Other negative responses observed in the herpetofaunal community can be addressed by changes to the macrohabitat (vegetation structure; Mushinsky, 1986), microhabitat (ground temperature, soil moisture, etc.; Blair, 1997; Knapp and Seastedt, 1986), or declines in prey abundances (most often due to lower abundances of arthropods in burned sites when compared to unburned sites; Evans et al., 1983). Annual burning is stressful to certain species, leading to local permanent emigration or local extinction. Past studies on effects of burning on herpetofauna have found variable responses among species. In tropical humid forests, herpetofaunal species richness was significantly higher in fire disturbed areas compared to undisturbed areas (Fredericksen and Fredericksen, 2002). Disturbances on small spatial scales at low frequencies may increase species richness by providing new habitat locally (McLeod and Gates, 1998; Mushinsky, 1985; Pilliod et al., 2003; Russell et al., 1999). However, there are also many studies which have reported that species respond negatively to burning (Cavitt, 2000; Rudolph and Burgdorf, 1997; Setser and Cavitt, 2003), or have a neutral response to disturbance by fire (Greenberg et al., 1994; Hannah and Smith, 1995; Woinarski et al., 1999). If immigration of new species or new individuals of an extirpated species is not possible, species richness levels decline.
Three caveats apply to the results of this study. First, we used an open population model based on the Jolly-Seber model. This model does not allow for heterogeneous capture probabilities across all different species, which may have been violated (Pledger et al., 2003). A violation of this assumption can lead to biases in parameter estimates, including underestimates of species richness (Huggins et al., 2003). A single sampling day per year was used to determine species occurrence, and some species may have been present but overlooked. Missing a species in multiple consecutive years may be due to permanent emigration or local extinction. Multiple survey days within a year would increase the probability of detection and allow for closed population modeling and estimation of species-specific capture probabilities. Last, our analysis was retrospective and an unplanned experiment with no replication. Lack of reference areas does not allow us to evaluate the possibility that declines in species richness were due to regional factors other than land management. Nevertheless, different decline rates in herpetofaunal species coincided with a switch in management practices (burning and stocking practices) from 1997–1999.
Management decisions based upon these results should be applied with caution, and future research with improved sampling methods is needed. Despite the caveats of our retrospective analysis, this 14-year data set remains one of the few long-term records of the herpetofaunal communities in the Flint Hills region of eastern Kansas.
The mark-recapture statistics used in this study provide an ideal method of analyzing community demographics of a multitude of ecosystems disturbed by a variety of human land-use practices.
Acknowledgments
Thanks to B. and J. Nitschke who own BJN Ranch and gave A. Volkmann permission to survey their land. Thanks to all other survey participants including N. Cubbage, J. Greider, R. Greider, J. Lent, E. McCarrier, J. Clark, R. Previtera, G. Volkmann, T. Volkmann, and S. Wiechmann. Financial support was provided by the Divison of Biology of Kansas State University.
